%% Estimate Structural VARs with Generalised Sign Restrictions
% by Jaromir Benes
%
% In this file, we show how to identify structural VARs based on sign
% restrictions, and other, more general, types of underdetermined
% restrictions.
%
% We estimate a simple reduced-form VAR model on made-up data. Defining a
% set of sign restrictions on various impulse responses to identify one
% particular shock, we randomly draw a given number of SVARs satisfying
% these restrictions. We then examine the impulse responses obtained from
% these SVARs.

%% Clear Workspace

clear;
close all;
clc;
irisrequired 20140602;
%#ok<*NOPTS>

%% Load Data
%
% Load the data made up in `make_up_data.m`. Remember -- these are not real
% economic series, just fake. So do not try to interpret the estimated
% correlations or impulse responses from an economic point of view.

d = dbload('VAR_Sign_Restrict_Data.csv');

d

%% Estimate Reduced-Form VAR
%
% Estimate a 2nd-order VAR using the three variables created above. Use the
% `estimate` function <?estimate?> with a database input and a list of
% variable names. The variable names, together with automatically created
% residual names, will be stored within the new VAR object `v` and used by
% other functions called below, such as `srf`. Note that the variable names
% and the residual names can be accessed through `get(v,'ylist')` and
% `get(v,'elist')`, respectively.

v = VAR({'R','PI','Y'});
[v,vdata] = estimate(v,d,Inf,'order=',2); %?estimate?

disp(v);

%% Discuss Identifying Restrictions
%
% Impose the following identifying restrictions (they all relate to the
% first shock; other shocks are not identified and will be ignored):
%
% * response of `R` (i.e. the first variable) to the first shock in the
% first period is positive;
% * response of `R` to the shock remains overall positive during the first
% two years (eight quarters), meaning that the sum of the interest rates
% over that period of time is greater than say 0.15;
% * response of `PI` (i.e. the second variables) to the first shock in the
% first two years is overall negative, meaning that the sum of inflation
% outcomes over that period of time is below -0.15;
% * response of `Y` (i.e. the third variables) to the first shock in the
% first two years is also overall negative.
%
% To write these restrictions formally, create an expression (as a text
% string) using references to individual impulse responses through variable
% named `S`. The `S(i,j,k)` element is the reponse of the `i`-th variable
% to the `j`-th shock in period `k`.

test_string = [ ...
    'S(1,1,1) > 0 ', ...
    '&& sum(S(1,1,1:8)) >= 0.15 ', ...
    '&& sum(S(2,1,1:8)) < -0.15 ', ...
    '&& sum(S(3,1,1:8)) < -0.15 '];

test_string

%% Identify Structural VAR
%
% Use the `SVAR` function with the option `'method='` set to
% `'householder'` to convert a reduced-form VAR to a structural VAR based
% on (generalised) sign restrictions passed; the `'householder'` refers to
% Householder transformations, see below. The restrictions themselves are
% passed in throught the option `'test='`.
%
% The procedure cycles over the following steps:
%
% # Randomly factorise the covariance matrix of reduced-form residuals to
% get a particular structural shock impact matrix (and hence a particular
% structural VAR);
% # Compute the impulse responses of that particular SVAR to the underlying
% structural shocks; these reponses are stored in an array named `S`.
% # Evaluate the test string for these impulse reponses. If the string
% evaluates to true the current structural VAR is kept. If it evaluates to
% false it is discarded.
% # Stop if the number of the accepted SVARs equals the number requested by
% the user (option `'ndraw='`) or if the number of cycles has already
% exceed the maximum number of iterations imposed by the user (option
% `'maxiter='`).
%
% The random factorisation used in IRIS is based on Householder
% transformations. The algorithm is described or used, for instance, in the
% following papers:
%
% * Rubio-Ramirez, J.F., D.Waggoner, T.Zha (2005). _Markov-Switching
% Structural Vector Autoregressions: Theory and Application_. FRB Atlanta
% 2005-27.
%
% * Fry, R. and A.Pagan (2007). _Some Issues in Using Sign Restrictions for
% Identifying Structural VARs_. National Centre for Econometric Research
% Working Paper 14/2007.
%
% * Berg, T.O. (2010). _Exploring the international transmission of U.S.
% stock price movements. Unpublished manuscript_.
% http://mpra.ub.uni-muenchen.de/23977.
%
% Comments on the code below:
%
% * <?sv?> `sv` is a structural VAR object with a number of alternative
% parameterisations.
%
% * <?B?> `B` is an array of the structural shock impact matrices that pass
% the test specified by the user in the option `'test='`.
%
% * <?count?> `count` is the number of iterations performed; here, do not
% restrict the maximum number of iterations (by setting `'maxiter'=Inf`),
% and hence it is the number of iterations needed to generate `N`
% successfully accepted structural VARs.

N = 500;

[sv,svdata,B,count] = SVAR(v,vdata, ...
    'method=','householder', ...
    'test=',test_string, ...
    'ndraw=',N, ...
    'maxIter=',Inf, ...
    'progress=',true);

sv %?sv?

size(B) %?B?

count %?count?

%% Plot Responses to Identified Shocks
%
% Simulate the impulse responses to the first shock (called `res_R` by
% default) in first four years. The option `'select='` limits the impulse
% response function `srf` only to that specific shock; otherwise all shocks
% would be simulated.

h = srf(sv,1:16,'select=','res_R','presample=',true);

h

[ff,aa] = dbplot(h,0:16,{'R','PI','Y'},'tight=',true);
ftitle('Responses to the identified shock');
ftitle('REMEMBER... these are just fake data, not real!', ...
    'location=','south');

%% Sort Parameterisations by Distance to Median
%
% Find the SVAR amongst all those identified in the previous step that
% somehow represents the typical impulse response. Use the SVAR function
% `sort`, which does the following:
%
% * Compute the user-selected impulse responses for each of the
% parameterisations within a SVAR object.
%
% * Find the median impulse reponses.
%
% * Compute the square distance to median for each parameterisation.
%
% * Re-order the parameterisations within the SVAR object so that the one
% closest to the median response comes first, etc. In other words, the
% first parameterisation will be "the most typical" one.
%
% Comments on the code below:
%
% * <?select?> Sort the SVARs by distance to median of the impulse reponses
% of all variable to the first shock in the first two years (eight
% quarters).
%
% * <?sv1?> The new object `sv1` is basically the same as `sv` except that
% the parameterisations are in different order now.

select = 'S(:,1,1:8)'; %?select?

select

[sv1,sv1data,index1,crit1] = sort(sv,svdata,select,'progress=',true); %?sortdata?

sv1 %?sv1?

figure();
plot(crit1);
grid('on');
title('(Sorted) Squared distance to median of response functions');
set(gca(),'xLim',[1,500]);
xlabel('Parameterisation #');

% ...
%
% Alternatively, if you're not interested in dealing with the SVAR data,
% call the `sort` function with an empty second input argument and waive
% the second output argument. Everything else, including the results,
% remains identical in <?sortdata?> and <?sortnodata?>.

[sv2,~,index2,crit2] = sort(sv,[],select,'progress=',true); %?sortnodata?

maxabs(crit1,crit2)

%% Plot Responses Closest to Median
%
% Compute and plot only the first `n=20` responses closest to median.
%
% Comments on the code below:
%
% * <?select?> Select the responses only to the first shock
% (named `res_R` by default); other shocks have not been identified.
%
% * <?link?> Link vertical axes in each of the graphs created here and
% those created above (if they still exist, i.e. have not been closed by
% the user in the meantime) so that the scale of y-axis is exactly the
% same.

n = 20;

h1 = srf(sv1(1:n),1:16,'select=','res_R','presample=',true); %?select?

h1

[ff1,aa1] = dbplot(h1,0:16,{'R','PI','Y'},'tight=',true);
ftitle(sprintf('The %g responses closest to median',n));
ftitle('REMEMBER... these are just fake data, not real!', ...
    'location=','south');

if all(ishandle(aa{1}))
    for i = 1 : length(aa{1})
        linkaxes([aa{1}(i),aa1{1}(i)],'y'); %?link?
    end
end

%% Help on IRIS Functions Used in This File
%
% Use either `help` to display help in the command window, or `idoc`
% to display help in an HTML browser window.
%
%    help data/dbload
%    help VAR/estimate
%    help SVAR/SVAR
%    help VAR/sort
%    help VAR/srf
%    help maxabs
